Simulation of Material Properties Below the Debye Temperature: A Path-Integral 

Molecular Dynamics Case Study of Quartz 
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Classical and path integral molecular dynamics (PIMD) simulations are used to study a and 
P quartz in a large range of temperatures at zero external stress. PIMD account for quantum 
fluctuations of atomic vibrations, which can modify material properties at temperatures below the 
Debye temperature. The difference between classical and quantum mechanical results for bond 
lengths, bond angles, elastic modulii, and some dynamical properties is calculated and comparison 
to experimental data is done. Only quantum mechanical simulations are able to reflect the correct 
thermomechanical properties below room temperature. It is discussed in how far classical and PIMD 
simulations can be helpful in constructing improved potential energy surfaces for silica. 



I. INTRODUCTION 

Quartz is one of the most abundant and best studied 
minerals on Earth; the structure and many other physi- 
cal properties are well understood!!!. Classical molecular 
dynamics (MD) simulations have been particularly suc- 
cessful in connecting atomic interactions between silicon 
and oxgen atoms in the condenspd phase with structural 



and elastic properties of quartzc Jland other silica poly- 
morphs at finite temperaturcsou'Lrtl. At zero tempera- 
ture, this conneetkin has also been done by mere first- 
principle studiesO til or .ab-inito methods incorporating 
bulk-system informationu. 

Some of the numerical approaches reach a nearly per- 
fect agreement with experiment, especially for structural 
properties. None of the calculations, however, include 
quantum effects of the ionic motion. Quantum effects 
lead to an equilibrium structure that is different from 
the "classical" equilibrium structure as soon as the in- 
teractions between the atoms constituting a crystal are 
not purely harmonic. Path integral simulations are a 
convenient, tool to compute such quantum mechanical 
effectst3~tJ, e.g., the zero-temperature lattice constant 
a(Ne cla ) of classical Neon is 0.18 A smaller than the 
one of Ife 22 , which again is about 0.19% smaller than 
a(Ne 20 ).ll 2 J Treating atomic motion quantum mechani- 
cally even effects ccwalent bonds p. systems such as crys- 
talline polyethylenell3 and siliconli-3. 

Quantum effects are certainly less strong in silica than 
in rare gas solids such as in neon. Thermal expansion 
can nevertheless be observed for a-quartz below the De- 
bye temperature and quantum effects have to be taken 
into account properly, if we want to relate simulations 
and experiments in a meaningful way. Correct quantum 
simulations will converge to a zero expansion coefficient 
as absolute zero is approached, while classical MD sim- 
ulations will (nearly) always result in a finite expansion 
coefficient even at T = K thus violating the third law 
of thermodynamics. It is obvious that structure and elas- 
tic constants are reflected well by the model potential if 
these properties have been used to fit the free parame- 



ters of the potential energy surface - as done in the case 
of the so-called BKS potential! It is therefore a much 
harder test for a model potential surface to yield the 
correct thermal expansion at low temperature, because 
the anharmonic interactions must be reflected accurately 
by the model potential. Overall, ab-initio calculations, 
MD simulations, and experimental data have partially 
achieved such good agreement that estimating quantum 
effects will play an important role in determining the real 
merits of a potential energy surface. 

In this study, path-integral molecular dynamics 
(PIMD) are used to determine the quantum effects on 
physical properties of quartz. Due to the large compu- 
tational demand of path integral simulations, we confine 
ourselves to the use of only one potential energy sur- 
face, namely the BKS potential! It has been particu- 
larly successful in reproducing silica properties not—only 
of a-quartzoa, but|-aJsa,of other silica polymorphscm and 
of the glassy stateli3~LD. It is therefore rather plausible 
that the BKS potential predicts accurately the shifts from 
classical results to quantum mechanical results. 

Recently, /3-quartz and /3-cristobalite have been inves- 
tigated by paeans of path integral Monte Carlo (PIMC) 
simulationstU. The PIMC study, however, was only done 
for phases that are stable at temperatures above 800 K. 
Quantum effects are relatively small at such high tem- 
peratures. Moreover, the resolution of our PIMD simu- 
lations exceeds that of the emploied PIMC algorithm by 
orders of magnitudes. E.g., the PIMD approach makes it 
possible to calculate the ground-state equilibrium lattice 
constants with a resolution of more than 0.001 A for a 
given potential energy surface, while the PIMC simula- 
tions had an uncertainty of typically 0.1 A. This improve- 
ment in resolution is possible despite a strong reduction 
in CPU time. Furthermore, arbitrary parallelepiped sim- 
ulation cells are permitted in this PIMD study allowing 
to calculate all elastic constants. The PIMC studies were 
confined to orthorhombic geometries and did not allow 
calculations of any elastic constants due to large statisti- 
cal error bars. 

The remainder of this paper is organized as follows: In 
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Sec. [IJ, the PIMD method used in this study is described 
along with some specific, technical details. In Sec. [II, re- 



sults are presented for structural data, elastic constants, 
and dynamic properties. Quantum mechanical results 
are compared to classical simulations and experiment. 
Conclusions are drawn in Sec. IV. 



II. METHOD 

A. Path Integral Molecular Dynamics in the 
Constant Stress Ensemble 

Although path integral Monte Carlo (PIMCVis usu- 
ally used to estimate quantum effects in. solidsc9~E3, path 
integral molecular dynamics (PIMD)c^l arise as a more 
natural choice for long-range Coulomb interactions. The 
Ewald sumEj has to be evaluated only once per time step 
in PIMD as opposed to each single local move in PIMC. 
Moreover, global moves in which the shape and size of the 
simulation cell are varied, are at basically no extra cost 
in molecular dynamics simulations, while Monte Carlo 
requires the evaluation of the net energy for each single 
trial move of the strain tensor. 

In the piath integral formulation of quantum statistical 
mechanics^, a quantum point particle at temperature T 
is represented by a closed polymer at temperature P T, in 
which adjacent beads interact via harmonic springs. The 
stiffness of the springs increases with decreasing thermal 
de Broglie wavelength A that a free particle would have 
at temperature P T. Studying this model in a molecular 
dynamics simulation, would require small time steps in 
the quantum limit P — > oo, if the dynamical masses of 
the polymer beads were chosen to be identical with the 
physical mass m of the quantum particle. However, it is 
possible to adjust all intra-molecular vibrations to similar 
time scales if the equations of motion are expressed in a 
convenient representation and appropriate "dynamical" 
masses are attributed to the beads.E3 In this study, the 
coordinates are represented in terms of eigenmodes of the 
free particle. The dynamical masses rh^, attributed to the 
motion of the eigenmode lu are usually chosen according 
to fh^jra = /ce/(&e + k u ), with the coupling of an 
atom to its lattice site in the Einstein model of solids 
and k u the spring constant associated with eigenmode 
lu. This choice of dynamical masses allows for efficient 
sampling of all degrees of freedom, because all modes 
move essentially on the same time scale. The dynamical 
mass of the center-of-mass motion of the polymer mo is 
of course identical with the real mass m. 

The motion of the simulation cell is constrained to 
symmetric strain tensors, but otherwise done as is in the 
classical Parrinello-Rahman methoccj. The dynamical 
mass W associated with the motion of the simulation 
cell geometry is again chosen such that a typical oscilla- 
tion time of the box is close to a typical oscillation time 
of a silicon or oxygen atom. The choice of W merely con- 



trols the efficiency of the sampling but leaves meaningful 
observables uneffected. 

One advantage of the Parrinello-Rahman method is the 
possibility|-tp determine all elastic constants at zero exter- 
nal stress.EH This is done by using appropriate relations 
between strain fluctuations and mechanical compliances. 
It is important to note that only isothermal strain fluc- 
tuations are accessible in PIMD simulations. A constant 
enthalpy simulation of the isomorphic classical represen- 
tation would not translate into conserved enthalpy of the 
quantum crystal. 

In principle, our method is closely related to a recently 
proposed PIMD scheme-, for constant-strain constant- 
temperature simulationsc3'E3. The special representation 
used here as well as omitting the th^Hnostat included 
in the equations of motions in Ref.c3'E3, anticipates to 
briefly review the final result. In order to do this, we 
represent the coordinate Ha of particle i at imaginary 
time t as a product of a scaled coordinate and the 
time-dependent (symmetrical) matrix h, which contains 
the shape and the volume V = det h of the simulation 
cell: 



Rita = h a @r 



(1) 



The values of ru a are constrained to values < ru a < 1 . 
The components of the metric tensor G are defined as 
G a f} = h afl hfj,f) where summation convention over Greek 
indices enumerating spatial dimensions is implicitly as- 
sumed. It is then convenient to express the equations of 
motion for the scaled coordinates in reciprocal Fourier 
space, namely in terms of coordinates 
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for which the motion of the free particles is diagonalized. 
Introducing k iw = Ami sin 2 (ttlu / P) / ({3h / P) 2 , allows to 
represent the equations of motion in a rather condensed 
form: 
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with Vij a two-particle interaction potential between par- 
ticle i and j and R.(y)t the vector connecting particle i 
and j at imaginary time t emploing minimum image con- 
ventions. 

Despite the weil-known disadvantages of Langevin- 
type thermostatsc-iE3, correlation times turned out to be 
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particulary small when all degrees of freedom (including 
the geometry of the simulation cell) were weakly coupled 
to a friction force-lipear in velocity and to a correspond- 
ing random force.E3 Chosing the damping term 7 of the 
Langevin dynamics to be 7 = 0.01 dt and dt = 0.04i c har 
with t c har the (smallest) characteristic time-scale of the 
system, systematic errors were made much smaller than 
statistical errors. Moreover, the ergodicity.-Broblcms, 
which is inherent to some PIMD algorithmscJ, can be 
most easily overcome with a Langevin thermostat. The 
average "dynamic" kinetic energy (T^ yn ) and the fluctu- 
ation of Xdyn in the PIMD approach described above, are 
utterly sensitive to bad choices of 7 and dt. In the cor- 
rect limit, one obtains (Td yn )/N = dksTP/2 per degree 
of freedom and the associated specific heat (fluctuation) 
(STl yn )/Nk h T 2 P 2 = dk B /2 with d the spatial dimension 
of the system. This sensitivity can be used to deter- 
mine the accuracy of the simulation resulting in reliable 
thermodynamic expectation values of other observables. 
Note that only independent components of the Fourier 
transform v iuJ need to be thermostated and considered 
for the calculation of the dynamic kinetic energy. 

B. Dynamical information 

Imaginary-time path-integral methods do not allow di- 
rect calculation of dynamical properties. While the com- 
plete dynamical information is obtained-jn imaginary- 
time correlation functions in principleEII, the inverse 
Laplace transform that one needs to carry out in order to 
assess real-time correlation functions is numerically un- 
stable. A generalization of the PIMD method, the cen- 
troid molecular dynamics (CMD) methodES, appears to 
give a more direct link between exact quantum dynamics 
and CMD. The basic idea of CMD is to propagate the 
center of mass of a PIMD quantum chain such that the 
internal degrees of freedom of the quantum chain are in 
their thermodynamical equilibrium. This makes the cen- 
ter of mass move on an effective classical potential sur- 
face. Time correlation functions can then be formulated 
in terms of centroid coordinates and averaged during the 
simulation. The Fourier transform of the centroid corre- 
lation function I c (u>) and the real quantum mechanical 
spectral function I(ui) are related linearly by 

I(u>) = n(Phw)I c (u) (5) 

with 

n(« = ^(l + coth^). (6) 

Eq. (||) is exact for harmonic systems as well as in the 
classical limit. Recently, much progress has been made 
in formulating rigorous relations between the dynam- 
ics of path-integral centroid variables and true dynam- 
ics as well as systematic correction to the CMD method 
sketched aboveS 



C. Model Specific Information 



In the present simulation, Vij ki Eqs. (|3|,[|) corresponds 
to the BKS interaction potentialu between particles i and 
j including the Coulomb energy among other contribu- 
tions. For the evaluation of the Ewald sum in arbitrar- 
ily shaped parallelepiped simulation cells, a recently pro- 
posed algorithmic was used and constrained to symmet- 
ric matrices h resulting in a 30% reduction of CPU time 
to evaluate the Ewald sum. The non-Coulombic interac- 
tions were cut off at a distance r c = 9.5 A. 

In the following, it will be distinguished between classi- 
cal and quantum mechanical results. Classical results are 
obtained by simply chosing P = 1 in a PIMD simulation. 
Exact quantum results require taking the limit P — > 00 in 
principle. It is well known that the so-called primitive de- 
composition of the density matrix, upon which the PIMD 
algorithm is based, invokes systematic errors in the parti- 
tion functioa and observables that vanish proportionally 
to 1/(TP)t13. It is therefore important to chose P large 
enough so that quantum effects are well reflected. On 
the other hand, P should not be too large, which would 
result in large statistical error bars. 

In order to assess at what values of TP one can expect 
convergence to the quantum limit, the average potential 
energy (V pot ) is calculated for various values of P at the 
temperature T = 300 K. The results are shown in Fig. |]. 
For PT > 1200 K, extrapolation to the quantum limit is 
possible with correction in the order 1/(TP) 2 . Combi- 
nations satisfying PT > 4, 800 K basically correspond to 
the quantum limit. 

Depending on the property of interest, it is necessary 
to extrapolate to the quantum limit with a 1/ P 2 correc- 
tions. This concerns particulary structural data, such as 
lattice constants, which can be obtained with high res- 
olution. Elastic constants and dynamical information, 
however, are plagued with relatively large statistical er- 
ror bars. For these observables, statistical uncertainties 
are much larger than systematic deviations and results 
obtained with PT > 2, 400 K are referred to as quantum 
limit. 

The total number of atoms used in the simulations 
was N = 1080. The linear box dimensions typically are 
24.9 A, 25.9 A, and 21.9 A along the a, b, and c axis, re- 
spectively. In order to obtain elastic constants with an 
accuracy of about 3 GPa, 60,000 MD steps of length 1 fs 
have to be performed, which takes about one day on an 
Intel II processor. For quantum simulations, the numer- 
ical effort has to be multiplied with Trotter number P. 



III. RESULTS 
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A. Structural Properties 



Two temperature regimes are particularly interesting 
to study, namely room temperature and temperatures 
near absolute zero. We first start with a discussion of 
quantum effects at room temperature, T = 300 K. This 
temperature is already well below the Debye temperature 
Td of quartz. Td of a-quartz as determined by specific 
heat measurements is a strongly temperature-dependent 
functionEEl: At T — 0, Td ~ 550 K, while at room temper- 
ature Td « 1, 000 K. Some quantum effects can therefore 
be expected to be relatively strong at room temperature. 
In Fig. ^, one can see that the main quantum effect at 
room temperature is the quantum mechanical freezing of 
the Si-0 bond. This freezing is reflected by the fact that 
the distribution psio( r ) of the Si-0 bond is much broader 
for the quantum mechanical study than for the classical 
study (Fig. ||a). In fact, the quantum mechanical psiO (?") 
barely changes its form when the temperature is lowered 
further. This is an indication that Si-0 bonds are in 
their quantum mechanical ground state. On the other 
hand, Si-O-Si bond angles as well as O-Si-0 bond angles 
do not differ considerably between classical and PIMD 
simulations (Fig. Bp). 

Fig. H and Fig.^| show the effect of the quantum me- 
chanical ionic motion on simple structural properties. 
The average Si-0 bond length rsio is shown in Fig. || 
as a function of temperature, while the lattice constants 
a and c are shown in Fig. [|. In all cases, it is noticeable 
that the quantum mechanical values are larger than the 
classical equilibrium lengths. The effects are relatively 
small, but clearly within the resolution of the simula- 
tions. While r Si0 only differ by 0.19% at T = 150 K, the 
lattice constants differ by 0.35% in the case of both the 
a axis and c axis of a-quartz. This means that the "ex- 
cess" quantum volume can be attributed to both the SiO 
bond length and_cwantum fluctuations of the so-called 
rigid unit modesOEZl. 



1 redirect com- 
^ The dif- 



In the case of the lattice constants, Fig. 
parison can be made to experimental dats 
ference in the lattice constants between quantum me- 
chanical calculations and experiment is about 0.06 A for 
the a axis and 0.07 A for the c axis. This difference is 
much larger than the discrepancy between classical and 
quantum mechanical simulations. It is, however, obvious 
that only quantum mechanical simulations reflect quali- 
tatively the right low-temperature behaviour: While clas- 
sical simulations result in a finite expansion coefficient 
even at absolute zero, the PIMD simulations lead to a 
vanishing expansion coefficient. 

There is, however, a very good agreement in the low- 
temperature thermal expansion along the c axis between 
PIMD simulation and experiment up to about 600 K, 
which shows that even anharmonic effects are well re- 
flected by the BKS potential as long as the system is 
still far away from the a-(3 transition. Near the transi- 
tion, however, the agreement becomes significantly less 



good. The "jump" in c at the transition is absent in 
the simulation. Thus, an important ingredient in the 
potential energy surface is missing. Note that classical 
simulations based on both the BKs and the TTAM po- 
tential do not reflect the anomaly fa the c/a ratio, which 
has been observed experimentallyO. Thus, while har- 
monic and low-temperature anharmonic effects are well 
described by the BKS potential, there seems to be the 
need to reflect many-body effects as well. 

Fig. H and Fig. || give detailed information on the bond 
angle distribution and their quantum effects. In Fig. |^, 
the Si-O-Si and O-Si-0 angle distributions are shown ex- 
emplarily at a temperature T = 80 K. Fig. || confirms 
the picture that the local structure in a-quartz does not 
correspond to tetrahedra. This can be concluded from 
the existence of two peaks in the classical O-Si-0 bond 
angle distribution and a broadened shoulder in the quan- 
tum mechanical O-Si-0 bond angle distribution. The dif- 
ference in quantum mechanical and classical mean bond 
angles is rather small, yet, noticeable, e.g., the classical 
average O-Si-0 bond angles approaches the ideal tetra- 
hedra angle of 109.471° much closer than the quantum 
mechanical simulation. The difference of about 0.04° be- 
tween the two approaches can be clearly resolved. The 
effect is much larger for the Si-O-Si bond angle, namely 
0.2° at T = 80 K, but less obvious (Fig. |b) because of 
the strong temperature dependence of (asiOSi)- 



B. Elastic properties 

Just like other properties can elastic constants be ex- 
pected to differ between classical and quantum mechan- 
ical treatments. In order to calculate classical constants 
at zero temperature, it is sufficient to calculate the second 
derivative of the ground state (potential) energy with re- 
spect to the stress tensor, resulting in the so-called Born 
expression for elastic constants]^ 



C aP =d 2 (V(T,e))/de a dep, 



(7) 



where the derivative is evaluated at zero strain. At finite 
temperatures, it is not sufficient to generalize this expres- 
sion by simpling taking the thermal expectation value p£ 
the right hand side. It has been pointed out corrective!! 
that the free energy surface 3-(T, e) should be considered 
instead of V(T, e). This generalization leads to differ- 
ent estimators of the elastic constants when evaluated 
in the (NVT) ensemble. The main effect of this gen- 
eralization is that fluctuations of the stress tensor need 
to be considered on top of the Born term described in 
Eq. 0. These fluctuations usually lead to a reduction of 
the elastic constants. Unlike classical fluctuation terms, 
quantum mechanically calculated terms will not vanish 
as the temperature approaches absolute zero. Among 
other effects, this will lead to different elastic -constants 
for quantum mechanical and classical systemscll. In the 
case of silicates, however, it turns out that it is more effi- 
cient to calculate elastic constants Cy by exploiting the 
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relational between dj and the thermal fluctuations of 
the strain tensor. 

Experimental, classical, and quantum mechanical elas- 
tic constants are compared in Fig. 0. Elastic constants 
can be expected to show larger (relative) quantum cpt- 
rections than lattice constants and heat of formationEJ. 
For quartz, the reduction of about 5 GPa in C33 seems 
to be the most dramatic effect. At T = 300 K, classical 
and quantum mechanical elastic constants agree within 
the statistical error bars. Below 300 K, the classical C33 
shows a stronger temperature dependence than the quan- 
tum mechanical C33. This effect should be taken into ac- 
count when trying to optimize potential energy surfaces: 
C33 predicted by the force field parameters for T = K 
should be a little larger than C33 measured at a temper- 
ature of 300 K. For other Cy the same comment applies 
in principle, but quantitatively, the effects are less dra- 
matic. 



C. Dynamical Properties 

As a generic dynamical property we consider the (clas- 
sical) inverse-mass weighted momentum autocorrelation 
function C(t) 



(8) 



where t denotes the real time. C(t)'s Fourier transform 
C(lu) can be used to define an effective density of states 

9<ss{v) 



Nk B Tn{(3hv) 



(9) 



with n(J3hv) being introduced in Eq. (ph. g c g(v) is iden- 
tical with the real density of states (DOS) if the har- 
monic approximation is valid. C(u>) and hence the effec- 
tive DOS can be exactly related to the imaginary-time 
correlation function G(t) 



Ri{T)-Ri(0) 
via the two-sided Laplace transform 

/oo 
dcu exp(-Tw0/2) 
-OO 



(10) 



cosh < Tilj 



cosh 



hup 



(ii) 



Note that the imaginary time r has to be considered 
within the interval < r < (3. Outside of this interval, 
imaginary-time correlation functions are repeated peri- 
odically. 

Eq. (|ll|) is useful to check the validity of the cen- 
troid molecular dynamics (CMD) method and hence to 



establish the validity of spectral functions as obtained 
by CMD. If C{t) is determined in terms of the mass- 
weighted autocorrelation function of the centroid veloc- 
ities and use is made of Eq. (||), the effective DOS and 
hence G(t) can be estimated in terms of centroid dy- 
namics. If CMD is applicable, G(r) as obtained by di- 
rect sampling and G(r) as estimated via CMD have to 
agree. As shown in Fig. |[ the agreement is perfect within 
our statistical error bars for a-quartz at very low tem- 
peratures. Of course, this agreement could be expected 
as the dynamics are dominated by the harmonic inter- 
actions unlike the thermal expansion coefficients. Note 
that a purely classical simulation leads to a similarly 
good agreement. CMD and classical velocity autocor- 
relation functions can barely be distinguished in the case 
of quartz. 

Fig. H shows the density of states as calculated via 
the centroid PIMD. Of course, a degree of (meaningful) 
complexity in a spectrum such as shown in Fig. can 
never be obtained by inverting imaginary-time correla- 
tion functions as shown in Fig. |[ Thus, centroid PIMD 
are a useful tool to obtain DOS of silica. 



IV. CONCLUSIONS 

This study shows that path integral molecular dy- 
namics (PIMD) are an efficient tool to calculate low- 
temperature properties of solids even if the complexity 
is larger than in rare gas crystals or other monoatomic 
solids. PIMD turns out to be particulary useful (as com- 
pared to path-integral Monte Carlo) when long-range 
forces have to be evaluated such as it is the case for the 
simulations of silica. Structural properties can be eval- 
uated with high resolution and the shift from properties 
that are obtained if atomic motion is treated classically 
to the "real" quantum mechanical properties can be as- 
sessed. This shift can also be calculated for elastic con- 
stants, which are notoriously hard to compute even in 
classical simulations. 

The result of the PIMD simulations anticipate that 
path integral techniques may not only become an impor- 
tant way of evaluating the merits and failures of potential 
energy surfaces, but PIMD might give valuable input to 
construct reliable model potentials. Here, the PIMD cal- 
culations of the thermomechanicaLproperties of a-quartz 
were based on the BKS potentials. The construction of 
the BKS potential was pioneering in the sense that ab- 
initio calculations were combined with bulk properties 
in order to fit the free model parameters. In the latter 
part, lattice constants and elastic constants were calcu- 
lated for a classical system at T = K from the (fit) 
parameters and adjusted such that agreement with ex- 
perimental "quantum mechanical" (finite temperatures) 
data was optimum. The PIMD results in combination 
with the classical MD results presented in this paper, 
show that this part of adjusting the parameters of the 



5 



BKS potential allows for further optimization. Of course, 
one can not necessarily expect to find a two-body poten- 
tial energy surface for silica that describes interactions 
much better than the BKS potential. 
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FIG. 2. a) Probability density p(r) to find an oxygen 
atom in a distance r from a silicon atom, b) O-Si-O (left) 
and Si-O-Si (right) bond angle distribution function p(a). 
Solid lines reflect classical simulations, dashed lines represent 
quantum mechanical simulations. Temperature T = 300 K. 
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Statistical error bars are about 2 GPa. 
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